from __future__ import print_function
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

start = time.perf_counter()


def compute_error_for_line_given_points(theta, x, y):
    theta = np.array(theta)
    y_predict = np.dot(x, theta)
    totalError = np.dot((y_predict - y).transpose(), (y_predict - y))
    return totalError / float(len(x))


# 梯度算法
def step_gradient(theta, x, y, learning_rate):
    gradient = (2 / float(len(x))) * np.dot(np.array(x.transpose()),
                                            np.array(np.array(y) - np.dot(np.array(x), np.array(theta))))
    # 我是使用y-y_predict，所以我的导数都是负数，上面算的是负偏导数的和，因此下面使用加法
    theta = theta + learning_rate * gradient
    return theta


# 开始迭代
def gradient_descent_runner(x, y, theta, learning_rate, num_iteration):
    for i in range(num_iteration):
        theta = step_gradient(theta, x, y, learning_rate)
    return theta


def run():
    x = pd.read_excel("E:\本科教学\大三下\数学建模\新疆大学数学建模校赛2021年度赛题及说明\A题\data.xlsx", header=None)
    x = np.array(x)
    y = np.array(x[:, 5]).reshape(len(x), 1)
    x = np.vstack((np.vstack((x[:, 1], x[:, 2])))).transpose()
    x = np.c_[np.ones(len(x)).reshape(len(x), 1), x]
    learning_rate = 0.00000005
    theta = np.zeros(len(x.transpose())).reshape(len(x.transpose()), 1)
    num_iterations = 100000
    print("Starting gradient descent at theta = {0}, error = {1}"
          .format(theta, compute_error_for_line_given_points(theta, x, y)))
    print("Runing")
    theta = gradient_descent_runner(x, y, theta, learning_rate, num_iterations)
    print(theta)
    print("After {0} iterations theta = {1}, error = {2}".
          format(num_iterations, theta, compute_error_for_line_given_points(theta, x, y)))
    return theta  # 返回theta的值画图


if __name__ == '__main__':
    x = pd.read_excel("E:\本科教学\大三下\数学建模\新疆大学数学建模校赛2021年度赛题及说明\A题\data.xlsx", header=None)
    x = np.array(x)
    y = np.array(x[:, 5]).reshape(len(x), 1)
    x = np.vstack((np.vstack((x[:, 1], x[:, 2])))).transpose()
    x_bias = np.c_[np.ones(len(x)).reshape(len(x), 1), x]
    theta = np.zeros(len(x_bias.transpose())).reshape(len(x_bias.transpose()), 1)

    z = x_bias.dot(theta)
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(x[:, 0], x[:, 1], y, s=10, c="black")
    xx = np.arange(0, 300, 0.5)
    yy = np.arange(0, 250, 0.5)
    X, Y = np.meshgrid(xx, yy)
    theta = run()
    Z = theta[0] + theta[1] * X + theta[2] * Y
    # 作图
    ax.plot_surface(X, Y, Z, cmap='rainbow')
    # 添加坐标轴(顺序是Z, Y, X)
    ax.set_zlabel('Z', fontdict={'size': 15, 'color': 'red'})
    ax.set_ylabel('Y', fontdict={'size': 15, 'color': 'red'})
    ax.set_xlabel('X', fontdict={'size': 15, 'color': 'red'})

    plt.show()

end = time.perf_counter()

print("运行耗时：", end - start)
